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ABSTRACT 


The analysis of transient signals using classical techniques is frequently not 
satisfactory. The Fourier analysis is based on stationarity of the signal, and transients are 
non-stationary. A new technique for the analysis of this type of signal, called the Wavelet 
Transform, was applied to artificial and real signals. A brief theoretical comparison 
between the Short Time Fourier Transform and the Wavelet Transform is introduced. A 
multiresolution analysis approach for implementing the transform was used. Computer 
code for the Discrete Wavelet Transform was implemented. Different types of wavelets to 
use as basis functions were evaluated. 
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I. INTRODUCTION 


It is believed that different transmitters have unique signatures, that is, the 
transient or turn on-tum off response is unique to a given transmitter or at least to a 
particular make. In this thesis we will examine this issue. The resulting process algorithms 
are also useful in other transient signal applications. We will make use of a new analysis 
technique called the Wavelet Transform. 

Analysis of stationary signals is a very well studied subject. The researcher can 
choose among a variety of methods and apply the best technique. However, the study of 
non-stationary signals, like transients and signals that show a behavior with sharp 
transitions, presents a challenge. The Fourier Transform method is not appropriate 
because it uses as basis function the complex exponential that extends over infinite time, 
while transients are of short duration. In 1946, Gabor [Ref. 1] proposed a modification to 
the Fourier Transform. A sliding time-window was introduced to the original transform 
equation with the intention to look into the signal by small portions. Assuming that this 
small portion has a stationary behavior, the analysis of the signal by the Short Time 
Fourier Transform (STFT) is appropriate. Other methods such as the Instantaneous Power 
Spectrum and the Wigner-Ville Distribution [Ref. 11], have also been developed but have 
similar limitations as the STFT. 

The purpose of this thesis is to investigate the use of the Wavelet Transform for 
the estimation / classification of certain transient signals. The results may be adapted to 
problems related to acoustic and electromagnetic transients. 
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II. TRANSIENT SIGNAL ANALYSIS 


A. THE SHORT TIME FOURIER TRANSFORM ( STFT) 

1. The Continuous STFT 

The Standard Fourier Transform of a signal is defined as 


00 

X(Q) = J jc (t)e-J Qt dt, 

-00 


( 2 . 1 ) 


or as 


= <*((), 


( 2 . 2 ) 


Equations 2.1 and 2.2 show that the Fourier coefficients are computed as inner products 
of the signal with a complex exponential basis function of infinite duration. It decomposes 
the signal into frequency components. If the signal components vary sharply in the time 
domain, after being transformed these local irregularities are spread out over the entire 
frequency domain loosing their localization. This approach may not be suitable if the signal 
of interest is non-stationary. However, even non-stationary signals can be assumed to have 
a stationary behavior for a short time. Based on this assumption, Gabor [Ref. 1 ] proposed 
the modification of the Standard Fourier Transform. A window g(t) of short duration and 

centered at location r was inserted in Equation 2.1 as follows 


X(C1) = J x(t)g(t-i)e-fi“dt. 

-00 


(2.3) 


Equation 2.3 is called Short Time Fourier Transform ( STFT ). If the analysis window 
g(t) is Gaussian the STFT is called the Gahor Transform. The analysis window g(t) is 
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considered to be part of the specification of the STFT. The short-time section of x(t) is the 
product x(t)g(t-r) and it becomes clear that changing the analysis window will, in general, 
change all the short-time sections and consequently the STFT [Ref. 3], A major drawback 
of the STFT is that g(t), as a fixed duration window, is accompanied by a fixed frequency 
resolution, allowing only a fixed time-frequency resolution [Ref. 2], This is known as the 
Heisenberg inequality that states that for a given transform pair 

(2.4) 

where er t and a a are the root mean square of the variance of the signal over time and 
frequency respectively. The variances are given by 

00 

J* t 2 \g(t)\ 2 dt 

= -. (2.5) 

J \g{t)\ 2 dt 

-00 

00 

J fl 2 lG(ft)l 2 d!Q 

= ^-. (2.6) 

J |G(Q)I 2 <£2 

-00 

Hence, either time or frequency resolution can be enhanced at the expense of a poor 
resolution in the other domain since there is a lower bound for the product given by 
Equation 2.4. Note that the time-frequency resolution is fixed for the whole 
time-frequency plane once a window has been selected. Figure 1 shows the tiling of the 
time-frequency plane for the STFT, where each rectangle represents the basic resolution 
cell, given by G/CTq . 
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Figure 1. Tiling of the Time-frequency Plane for the STFT 
2. The Discrete-Time STFT 

The discrete-time version of the STFT can be defined as 

X(k,n) = 2 g(n - m)x{m)e~ J ^ km . {2.1) 

m=—oo 

The analysis window g(m) is considered to be part of the specification of the STFT, and 
its selection dictates the trade-off between frequency and time resolution. The 
discrete-time STFT can be viewed as a set of Discrete Fourier Transforms ( DFT ) of 
short duration. It can also be viewed as the output of a filter bank. 

Rewriting Equation 2.7 we have 

X(k,n) = £ [x(m)e-^]g ( w - m ). (2 8 ) 

m=-oo 
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- ~b? 

Equation 2.8 is in the form of the convolution of the sequence \x(m)e 1N ], with the 
sequence g(m). We can then rewrite Equation 2.8 as 

X(k, n) = [x{n)e~ i ^ kn ] * g(n), (2.9 ) 

where the symbol " * " denotes the convolution operator. Thus, the discrete-time STFT 
can be viewed as a collection of sequences, each corresponding to the frequency 
components of x(n) falling within a particular frequency band. Figure 2 illustrates a filter 
bank operation, where each filter acts as a bandpass filter centered at a selected frequency 
[Ref. 3], 


x(n) 



Figure 2. The Discrete-time STFT as a Filter Bank 
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To overcome the time-frequency resolution constraint of the STFT a new 
technique was devised to create a basis function to allow a multiresolution analysis. This 
new technique is known as the Wavelet Transform. 


B. THE WAVELET TRANSFORM 

To obtain a multiresolution analysis we must allow both resolutions a t and cr n to 
vary in the time-frequency plane. When the analysis is viewed as a filter bank, the time 
resolution <J t increases with the central frequency of the analysis filters, consequently the 
frequency resolution a a is inversely proportional to £2, or 


Q “ 


c. 


( 2 . 10 ) 


where c is a constant. The analysis filter bank is then composed of bandpass filters with 
proportional bandwidth (the constant Q analysis). Figure 3 shows a division of the 
frequency domain for the STFT and the Wavelet Transform. 



i_i_i_i_i_ 

7t/ 8 tc/ 4 2 n © 

Figure 3. Division of the Frequency-domain for the STFT and the Wavelet Transform 
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The Wavelet transform is founded on dilations and translations of a basis function 
known as the mother wavelet. The family derived from this basis function is of the form 

= (2ii) 

where a is the scaling parameter and b is the shifting or translation parameter. If a is large 
then Wab(t) becomes stretched and corresponds to a low frequency function and if a is 
small yjt) becomes compressed and corresponds to a high frequency function. The 
dilations and translations present in the wavelet function enable the narrowing and 
widening of the time-frequency plane. Thus, it is possible to adjust to the characteristics of 
the analyzed signal [Ref. 5]. Figure 4 shows the tiling of the time-frequency plane for the 
Wavelet Transform. Figure 5a shows a mother wavelet function while Figures 5b, 5c and 
5d show examples of its dilations and translations. 


frequency 



scales time 


Figure 4. Tiling of the Time-frequency Plane 
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Figure 5c. a = 2, b = 1 Figure 5d. a = 1, b = 1. 

Figure 5. Haar Wavelet and Part of its Family 


1. The Continuous Wavelet Transform 

The Continuous Wavelet Transform (CWT) can be defined as 


C.WT x (a,b) = -ji 

—00 


( 2 . 12 ) 


where is given by Equation 2.11 and the term \f Ja is the wavelet energy 

normalization factor. Equation 2.12 can be written as 
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CWT x (a , b) = 


(2.13) 


The inverse of the Wavelet Transform is given by 


00 00 

*(0 = £ / \cWT x (a,byv*(tff • 
-00 0 


where C v is given by 


00 



0 


(2.14) 


(2.15) 


and 'F(Q) is the Fourier Transform of yft). Equation 2.14 holds only if the wavelet 

function ) satisfies the admissibility condition. To be an admissible wavelet function, 
the term C y < oo, therefore y/^ (t) has zero mean or its DC value is zero ( IfO) = 0 ). The 

function y/^ (t) has compact support, hence it acts as the impulse response of a bandpass 
filter that decays very fast to zero [Ref. 2], 


2. The Discrete Wavelet Transform 

The parameters a and b of the mother wavelet yjt) are continuous variables 
making the Wavelet Transform very redundant and impractical. Both parameters can be 
discretized by sampling to obtain a set of discrete wavelet functions. However, if the 
sampling rate is too high, the transform is close to the continuous transform case and 
hence it might be redundant. If it is too low, it may not contain the necessary amount of 
information for reconstructing the original signal. The solution to the problem may be 
addressed by the use of frames. The study of frames was introduced by Duffin and 
Schaeffer [Ref. 6] and was also addressed by Daubechies [Ref. 7] Their work provided a 
general theory that sets the balance for the extreme cases of over and undersampling the 
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parameters of the Short Time Fourier Transform and the Discrete Wavelet Transform and 
proved the existence of numerically stable inversion procedures for both transforms. 


Let the sampling lattice be 

a = a~Q, (2.16) 

b = kboa~Q, (2.17) 

where j, k € Z . Equation 2.11 becomes 

Vab(t) = -j=- \|/(V) = (2.18 ) 

Substituting 2.16 and 2.17 into 2.18 we obtain 

M= 4 V(d 0 t - kb 0 ). (2.19) 


The set of wavelet functions defined by 2.19 can be considered as a discretized version of 
the continuous case and we can represent any signal or function jc(t) e L 2 (R) as 

x(t) = E S djJcVjk(t), j , k e Z, ( 2.20 ) 

j k 

where djjc is defined as 
00 

djjt = —■ J x{t)\y(d 0 t - kb 0 )dt. ( 2.21) 

-00 

The wavelet basis function cannot be considered a basis for L\R)\ it is called a frame. The 
frame does not satisfy Parseval's power theorem and the wavelet expansion may lead to 
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more than one solution. Daubechies [Ref. 7] showed that, the recovery of xft) is possible 
only if 


J k 


j,keZ, 


( 2.22 ) 


uu 

where A and B are the frame bounds for the frame 1} fjkif) and |M| 2 = J \x(t)\ 2 dt. The set 

—00 

defined by 2.19 is considered a frame only if vp(/) satisfies the admissibility condition. 
Then the frame bounds are constrained by 


OD 

A<-^— f 

bologon J 


mn)l 2 jn 


bologoo J |£2| 2 

-00 


<B, 0<A <B<co. 


(2.23 ) 


If Vj \Jjk(t) is inadmissible it leads to a diverging upper bound ( B = oo ). The inequalities 
hold for any a 0 and b 0 . Some definitions are needed before we proceed : 

> If A = B = 1 then the frame is named a tight frame; 

> If the removal of a function leads to an incomplete set, it is an exact frame; 

> An exact and tight frame constitutes an orthonormal basis [Ref. 7], 

For the set of functions defined by 2.19 to be considered a basis for L 2 (R) it has 
to satisfy the definitions above and, therefore, the tightest set vj/y*(/) is the orthonormal 

set. 

Let the set lbe a set of orthonormal wavelets then 


J y,u(t)yfk'(t)dt = | q 


J = f. * = ^ 
otherwise. 


(2.24 ) 
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Both indices are orthonormal in time within any scale j and also across all other scales 
[Ref. 2], 


For practical purposes a 0 = 2 and b 0 = 1, leading to what is called a dyadic or 
octave lattice [Ref. 2, 4, 7], Figure 6 shows the dyadic or octave sampling lattice. 
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Figure 6. Sampling Dyadic Lattice for the DWT 


If we substitute the values a 0 = 2 and b 0 = 1 in Equation 2.19 we have 


= 2 jf2 \y(2 J t - k ). 


(2.25) 


Mathematical and practical perceptions of the wavelet transform can be better understood 
by using the concept of multiresolution to define the effects of changing scales [Ref. 8], 
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3. Multiresolution Signal Analysis 

Multiresolution analysis of L 2 (R) is defined as a sequence of closed subspace V } of 
L 2 (R),j e Z that satisfy the following properties [Ref. 9]: 

> Containment 

...F_ 2 cF-icFocFicF 2 ... ; 

<- coarser finer -» 

> Completeness 

n Vj = {0} for ;'eZ, (i.e., there is no intersection of subspaces), 
U Vj = L 2 (R) for jeZ, (i.e., the union of subspaces is equal to L 2 (R))\ 


> Scaling 

x(t) e Vj ox(2 1) e Vj+ 1 for any function x e L 2 {R) \ 

> Shift or translation 

x(t) € Fo <=> x(t + k) e Vo for any function Jt € L 2 (R),k € Z; 

>■ A function cp(t) E Fo exists such that, for V/w € Z, the set defined by 

(Pjk(t) = 2 jn (p(2 J t — k) forms an orthonormal basis for Vj. (2.26) 

The functions defined by 2.26 are called scaling functions, since they are scaled versions 
of functions in L 2 (R) . From the above definition for multiresolution analysis, it can be said 
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that the function (p(7) in L 2 (R ) is the limit of the approximations cp/(7) e Vj as j 
approaches infinity [Ref. 5] or 


(p(0 = lim q>jk(t). 

j —>oo 


The variable j is called a scale factor. For j > 0, cp/fe(0 narrows, and consequently Vj is 
wider, being able to represent finer detail. If j < 0, cp jiff) stretches, Vj is narrower and 
the resolution is coarser. Let cp(V) c Vo be a scaling function; then its translations 
cp(f - k) span Vo . Let cp(2 1 - k) be dilations of cp(? -k) in Vo. We know that Vo a V\ 
then cp(2 t-k) span V\ . Therefore, (p(0 can be presented as a linear combination of 
translations of cp(2/) as 

cp(0 = Yj ho{k)(p{2t-k), keZ. ( 2.27) 

k 

Equation 2.27 is known as dilation or fundamental equation and ho(k) are the scaling 

function coefficients. We will see later that /io(£)can be interpreted as a FIR filter. The 
space L 7 (R) can be represented via a Venn diagram as shown in Figure 7. 
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Another way of graphically representing L 2 (R) is shown in Figure 8. 



L 2 (R) 


Note in Figure 8 that the areas denoted by Wo, W\, W2, W3 etc., are the difference 
spaces between (Fj, Fo),(F2, Fi), (F3, p2),(^4, F3), respectively. The idea can be 
extended to an arbitrary Wj. The spaces Wj are defined to be the orthogonal complement 

of the spaces Vj, with respect to V J+ \ [Ref 5 ] that is 

Vj+i = Vj © Wj, ( 2 . 28 ) 

Vj 1 Wj. ( 2 . 29 ) 

Where the symbol " © * stands for direct sum. This means that each element of V J+ \ can 
be written as the sum of an element of Wj and an element of F, . A subspace Wj 
provides the necessary detail information to go from subspace V } to V J+ \ . The direct sum 

of all the possibly infinite spaces is equal to L 2 (R) [ Ref. 10 ]. 
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Let lj f(t), as defined by Equation 2.25, be an element of any given subspace Wj. 
The subspace Wj is itself a subspace of Vj+\. Therefore, V| f(t) can be written as a linear 
combination of translations of (p(t) 

VJ/O) = Z h 1 (k)q(2t -k\ keZ, (2.30) 

k 

where h\(k) represents the wavelet coefficients that can be used as a high pass Finite 
Impulse Response (FIR) filter. These coefficients are related to the scaling coefficients, as 
they span two orthogonal spaces, they must also be orthogonal. Note that both equations 
are derived from a common term and differ only by the coefficients ho(k) and h\(k). 

Hence, orthogonality can be achieved by requiring that 

(ho(k),hi(k)) = 0 . (2.31) 

Assuming that the number of coefficients defining ho(k) and h \ (k) is even and 

hi(k) = (-\) k h 0 (N-l-k), (2.32) 

where N is the number of coefficients, then the inner product of the coefficients is equal to 


S ho(k)h\(k) = 2 h 0 (k)(-\) k h 0 (N- 1 -k). ( 2.33 ) 

b=0 k =0 

The expansion of 2.33 is equal to zero. Therefore, the space L 2 (R) is spanned by a linear 
combination of functions in the subspace Vo and functions in subspaces fV Q to IF . As a 
consequence, any function x(t) in L 2 (R) can be written as 
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x(t)= £ a k cp*(0 + S Z djk\\j jk (t). 

k =-oo y=l k-~<xj 


(2.34) 


Equation 2.34 is an expansion in terms of scaling and wavelet functions. The coefficients 
Clk and djk are the Discrete Wavelet Transform coefficients. They completely describe the 

original signal and are defined as: 


ak = (x(t), c p k (k)) = J x(t)(p k (t)dt. 


(2.35) 


□u 

djk = (x(tf \\tj k (t)) = f x(t)\\fj k (t)dt. 


(2.36) 


Where (p*(^) and vj \ijk(f) are real functions. If the scaling and the wavelet functions form 
an orthonormal basis, Parseval's power theorem holds. For the expansion defined by 2.34, 
Parseval's theorem is 


\\x{t)\ 2 dt- f) \a k \ 2 + £ f) \djk\ 2 . 

k=-€D f= ] k=- oo 


(2.37) 


If the functions are orthonormal and have compact support ( finite duration ), the desired 
time localization in all scales is achieved [Ref. 8, 10], 


a. The Scaling Function and the Scaling Coefficients 
The coefficients h Q (k) in Equation 2.27 have to be chosen carefully in 
order to generate scaling functions that have compact support (finite duration). They must 
satisfy some necessary conditions to ensure this compact support [Ref. 10, 11, 12], The 
development of these conditions are provided in Appendix A. 
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1. The sum of the scaling coefficients must be equal to two 


N-l 

2 


n =0 


h 0 (n) = 2. 


(2.38) 


2 . If a solution to 2.27 exists, and J (p(t)dt * 0, and J \(${t)\ 2 dt = 1 , and an 
integer number of translations of Cp(?) are orthogonal to each other, as defined by 

<(p0)(p*(0> = J <p(0q>(? - k)dt = 8(k) = 


f l,k = 0 

( 0, else 


then 


N~\ 

2 


n=0 


h 0 (n)h 0 (n - 2k) = 28{k) . 


(2.39) 


3. The sum of the squares of the coefficients is equal to two. 


Let k — 0, Equation 2.39 reduces to 


2 \ho(n)\ 2 = 2. 

n=0 


(2.40) 


4. The individual sums of the odd and even terms of h 0 (n) is equal to one 


S i +^2 = 2, 


where 


Si = I h 0 (2n), 

and 


N-\ 

S 2 - 2 ho(2n+\). 

n =0 


(2.41) 

(2.42) 


(2.43) 
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b. Scaling and Wavelet Functions Regularity 

The N number of scaling function coefficients has to be even. These 
coefficients have to satisfy the linear constraint determined by Equation 2.38 and the y 
bilinear constraints determined by Equation 2.39. This leaves y- 1 degrees of freedom 
for choosing them. That will guarantee the existence of the scaling function, and the set 
derived from it will be orthonormal [Ref. 10], 

The degrees of freedom in the scaling function coefficients determine how 
regular the function will be. Regularity or smoothness of a function is defined by the 
number of times that a function can be differentiated. In other words, a regular scaling or 
wavelet function can be defined as a function with compact support in time and reasonably 
localized in frequency. 

c. The Frequency Domain 

The Discrete Fourier transform of Equation 2.27 is defined as 


//(©)= £ h Q {k)e~^ k . 

k=~OD 


(2.44) 


It turns out that a solution to Equation 2.27 exists only if H(0) = 2, or the DC gain of the 
FIR filter equals two. 

If the integer translates of Equation 2.27 are orthogonal, then the 
constraint defined by Equation 2.39 is true, if and only if 

|//(io)| 2 + |//((0 +Jt)| J = 4. (245) 

Equation 2.45 is equivalent to a Finite Impulse Response (FIR) filter and is called a 
Quadrature Mirror Filter (QMF) [Ref. 10 , 18], This allows the implementation of the 
Discrete Wavelet Transform via a Filter Bank technique. 
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4. Filter Banks and the Discrete Wavelet Transform 

The multiresolution analysis introduced in the previous section can be implemented 
by applying a technique called Pyramid decomposition, or Mallat's algorithm [Ref. 13], 
The algorithm does not include the use of scaling or wavelet functions; only the scaling 
and the wavelet coefficients need to be considered. We derive the basic relationship 
between the expansion coefficients at a lower scale in terms of those at a higher scale. 
Reproducing the fundamental equation 

9(0 = £ /2o(&)(p(2f -k), JceZ. (2.46) 

k 

Next, we dilate and translate Equation 2.46 to obtain 

cp(2 j t - ri) = 2 ho(k)q>(2(2 J t -ri)-k) = 

k 

=£ h 0 (k)(p(2 J+] t - 2n - k) . (2.47) 

k 

Let m — 2n + k or k = m - 2n, substituting into Equation 2.47 

cp(2 H -ri) = Yjho(m- 2«)(p(2 y+1 / - m). (2.48) 

m 

Let Vj be the subspace spanned by (p(2 J t - k) then 

x(t) e Vj +1 => x(t) = £ a jk cp(2> +l t-k). (2.49) 

k 

Consequently x(t) can be expressed in terms of scaling functions and no wavelets. At one 
immediately lower scale, the wavelets will provide the necessary detail not available in 





that scale space and for a lower resolution the detail has to be accounted for. Therefore, 
the wavelets have to be included in Equation 2.49, which is now rewritten as 

*(0 = 2 a jk 2 jl2 y(2 J t - k) + 2 djk2 J,2 \\f(2 J t -k). (2.50) 

k k 

Where the term 2 jl2 is the normalization for the various scales, and the terms Cijk and dj k 
are the discrete wavelet transform coefficients. Note that x(t) belongs to Vj+\ and 
(p (2 J t — k) and \|/(2 J t — k) belong to Vj. If we substitute Equation 2.26 into Equation 
2.50 we obtain 

x(t) = 2 a jk (?jk(t) + 2 djkVjkit), k g Z . (2.51) 

k k 

The discrete wavelet transform coefficients can be written as 


00 

Ojk = J x(t)(p jk (t)dt. 
—00 


(2.52) 


Substituting Equations 2.26 and 2.48 into Equation 2.52 leads to 


GO 

a jk = | x(/)2 ho(m - 2A:)2 (/+1)/2 (p(2 7+1 t - m)dt. 


or 


00 



(2.53) 


Note that the integrand of Equation 2.53 is the discrete wavelet transform coefficient of 
the scale j + 1, which leads to the following relation 
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(2.54) 


cijk = 'Lho(m-2k)aj + \ yk . 

m 

The analogous relationship for djk is of the form 

djk = 'Ehi(m- 2 k)a j+ ijc . (2.55) 

m 

Therefore, the scaling and wavelet coefficients at a given scale j are obtained by 
convolving the coefficients at that scale with hoi-m) and and then decimating to 

produce the expansion coefficients at scale j-1 . This is equivalent of filtering the j th scale 
coefficients with two FIR filters. The impulse responses are ho(-m ) and hi(-m). After 
decimation or down sampling, the next coarser scaling coefficients are obtained [Ref. 10], 

Figure 9 shows the implementation of Equations 2.52 and 2.53 for three levels of 
analysis. The notation LP stands for the lowpass FIR filter (i.e., the weights ho(-m ) ) and 

HP stands for the highpass FIR filter (i.e., the weights h\(-m) ). 



Figure 9. Three levels of Wavelet Transform Analysis 
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The implementation of the algorithm was done using matrices. An example of a 
DWT using a four coefficient wavelet is given as follows. 

Defining a wavelet coefficient matrix A as 


h 0 (0) h 0 (l) h 0 (2) h 0 (3) 
h 1 (0)h 1 (l) h 1 (2)h 1 (3) 

h 0 (0)h 0 (l) h 0 (2)h 0 (3) 
h 1 (0)h 1 (l) h 1 (2)h 1 (3) 


h 0 (2) h 0 (3) 
h 1 (0)h 1 (l) 


h 0 (0)h 0 (l) h 0 (2)h 0 (3) 
h,(0) hj(1) hj(2) hj(3) 

ho(0)h 0 (l) 

h i (2) h,(3) 


The matrix A is a sparse matrix composed of scaling coefficients ho(m), and wavelet 
coefficients h \ (m) as non-zero entries. The blanks indicate zero entries. If we multiply a 
given signal by the matrix A, the first row generates one component of the signal 
convolved with the scaling coefficients. The following odd numbered rows perform the 
same operation. The even numbered rows generate the convolution of the signal with the 
wavelet coefficients. Note that there is an offset by two every two rows. This offset by 
two is the equivalent of a decimation operation. The last two rows wrap around like 
convolutions with periodic boundary conditions [Ref. 15], 

The DWT consists of applying the matrix A repeatedly to a signal with 
starting length N until a stopping criterion is reached. In this thesis, when the number of 
scaling coefficients was equal to one, we stopped. 

As an illustration, let x(t) be an arbitrary signal of length 16. Applying a 
DWT to it, we obtain: 
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Figure 10. DWT of an Arbitrary Signal x(t) 

Note that at each level matrix A is applied to half of the values of the previous iterations. 
Once the d coefficients are generated, they just propagate through all subsequent levels. 
The DWT can operate only at signals with length equal to a power of two (2 ",n e Z). 
These matrices operations are generally carried until there is only one DWT scaling 
coefficient a,,. 

The Matlab code used for implementing Mallat's algoritmh in this thesis plots the 
scales of the DWT, with a number of points equal to half of the total length of the signal 
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III. WAVELET SIGNAL PROCESSING 

A. INTRODUCTION 

The Discrete Wavelet Transform was implemented using the Mallat's Pyramid 
algorithm described in Chapter D.4. All programs for the DWT were written in Matlab® 
[Ref. 17] and are given in Appendix B. The processing phase was divided in two parts. A 
simulation phase where the DWT was applied to four artificial communication signals, 
and a testing phase where the DWT was applied to three types of real signals. The real 
signals were composed of transient signals created by metallic objects, radar pulses, and 
turn-on / turn-off transients of push-to-talk radios. Six types of mother wavelet functions 
were used to process the signals. 

B. MOTHER WAVELETS 

The Haar wavelet is the simplest type of wavelet basis function. It is very well 
localized in time, but the frequency localization is poor due to the discontinuities in the 
time domain. Figure 11 shows the Haar wavelet. 



0 0.5 1 1.5 2 

Figure 11. The Haar Wavelet 


Figure 12 shows the one-sided spectrum of the Haar wavelet and its dilations. 
Note that the spectral sidelobes produce a poor frequency localization. 
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0 0.5 1 1.5 2 2.5 3 

(radians) 

Figure 12. One Sided Spectrum of the Haar Wavelet and its Dilations 

The Daubechies' scaling and wavelet functions are the most popular functions for 
use in DWT. Daubechies' scaling coefficients satisfy all necessary conditions. The wavelets 
are highly regular. They have compact support in time and their use in the DWT leads to a 
reasonable localization in frequency. The regularity of this type of wavelets increases 
linearly with the length of the FIR filter [Ref. 2], Daubechies' wavelets are identified by the 
number of their coefficients, e.g., a four coefficient wavelet is called Fourth Order 
Daubechies or Daubechies four. In this thesis, the short notation Daub(n) is used to 
denote the n* order Daubechies filter. Figures 13, 14, and 15 show the fourth, twelfth and 
twentieth order Daubechies' wavelets, respectively. 



Figure 13. The Fourth Order Daubechies Wavelet 
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Figure 14. The Twelfth Order Daubechies Wavelet 



Figure 15. The Twentieth Order Daubechies Wavelet 


Note that as the order of the wavelet is increased it becomes smoother. However, if 
computation speed is a factor then a lower order Daubechies' wavelet might be preferred. 

Figures 16, 17, and 18 show the frequency response of Daub4. Daub 12 and 
Daub20 wavelet functions, respectively. Note in Figure 16 the sidelobes around 2.5 
radians. Observe that in the subsequent figures these sidelobes are gradually becoming 
smaller, as the order of the wavelet increases. 


29 








figure 16. One Sided Spectrum of the Daub4 Wavelet and its Dilations 



0 0.5 . 1 1.5 2 2.5 


radians 

Figure 17. One Sided Spectrum of the Daub 12 Wavelet and its Dilation; 



Figure 18. One Sided Spectrum of the Daub20 Wavelet and its Dilations 







The Quadratic 5-spline wavelets are piecewise quadratic functions. These wavelets 

are called the Battle-Lemarie bases [Ref. 14] and are referred to in this thesis as spline 
functions. The initial scaling function cp(Y) is compactly supported. However, the 
functions 2 J,1 q>(2 J t — k) derived from it are not orthogonal. The spline scaling coefficients 
do not satisfy all the constraints described in Subsection D.B.3.a. The use of the spline 
function in the DWT leads to a redundancy in the scales. Transient signal analysis may, 
sometimes, benefit from this redundancy. These functions have an arbitrarily high number 
of derivatives [Ref. 7], which makes them extremely regular. Figure 19 shows the four 
coefficients spline function, and Figure 20 shows the one-sided spectrum of its dilations. 



Figure 19. The Fourth Order Spline Wavelet 



Figure 20. One Sided Spectrum of the Spline Wavelet and its Dilations 
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The ad hoc wavelet is a trial-and-error wavelet function. The of purpose of using it 
was to find a matching wavelet function. This function would be matched to the signals 
characteristics and, therefore, enhances the detection of the signal when wavelet 
transformed. 




0 0.5 1 1.5 2 2.5 3 


radians 

Figure 22. One Sided Spectrum of the ad hoc Wavelet and its Dilations 

C. THE SIMULATION PHASE 

In this part, four artificial communication signals created by Yayci [Ref. 16] are 

used as test data. These signals have sharp transitions in time that are exploited via the 
wavelet transform. They are: Binary Phase Shift Keying (BPSK), On-Off Keying (OOK), 
Frequency Shift Keying (FSK) and Quadrature Phase Shift Keying (QPSK). They are 
shown in Figures 23 through 26. 
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Figure 26. QPSK Signal 

1. Application of the DWT to Artificial Signals 

The DWT was applied to each of the communication signals. The first experiment 
used noise free signals. Gaussian White Noise (GWN) was added to the signals with six 
levels of Signal-to-Noise-Ratio (SNR). The SNR levels were 20 dB, 15 dB, 9 dB, 6 dB, 3 
dB, and 0 dB. Figures 27 to 33 present the decomposition of the noise free signals. The 
order of the scales in these plots are reversed, i.e., the first DWT level appears at the end. 
Only the last four scales are shown because the other ones have few coefficients and little 
energy. The decompositions for the other levels of SNR are shown in Appendix B. 

The transition points of the signals are easily seen in the noise-free case but 
become harder to identify as the SNR decreases. Preprocessing the noisy signals using 
differentiation, Hilbert transforms and a combination of these two techniques did not seem 
to improve the identification of the transients. Differentiating a signal in additive GWN 
causes the scales to display false peaks. The increase in the order of Daubechies' wavelets 
improves the probability of identification of the transitions. The ad hoc wavelet performs 
equally well for all cases of GWN. The transitions on the BPSK and the OOK signals can 
be identified at lower SNR relative to the FSK and QPSK signals. This is, in part, due to 
the presence of sharper transitions in time. 
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2. Evaluation of the DWT for the Artificial Signals 

The ability of identifying the transients of the artificial signals embedded in noise 
was the performance measure for the wavelet functions. The true transient locations were 
known a priori. Therefore, a simple comparison between the peaks shown in the 
transform domain and the known transition instances was sufficient to determine how well 
the wavelet function performed. The figures used for comparison are shown in Appendix 
B. The ad hoc wavelet function clearly outperformed the other ones. This wavelet was 
chosen by trial-and-error, and it produced the strongest peaks at the transition instances. It 
is possible that its efficiency will be degraded when working with other types of signals. 
The spline wavelet achieved a reasonable result down to 6 dB of SNR. The twelfth order 
Daubechies’ wavelet presented acceptable results for all signals of interest. Hence, it was 
chosen to be the basis function for the analysis of the real signals during the test phase. 

D. THE TEST PHASE 

This section presents a description of the signals used and their respective 
transforms. Only the last four scales of the DWT are shown. The signals were divided into 
three groups of transients: sounds produced by metallic objects, radar pulses, and tum-on / 
turn-off transients produced by push-to-talk radios. These transient signals are embedded 
in noise with unknown power levels. 

1. Metallic Sound Transients 

The first group of transient signals was obtained from the sound files directory of 
the Digital Signal Processing group of the Electrical and Computer Engineering 
Department Computer Center. There is no documentation on how they were generated or 
recorded. The length of the transients was cut to the immediate lower power of two. 
Figures 34.a to 34.c show the time representation of sounds produced by a dropped 
metallic object on a metallic surface, a gong, and a latch being opened, respectively. 
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Figure 35 shows the DWT of the last four scales of the metallic sound transients using a 
Daub 12 wavelet. 

Each signal has a particular sound. The time representation of the gong and the 
dropped metallic object show that both signals have an exponential decaying. The sound 
produced by the latch is clearly different from the other ones, with no particular time 
behavior. 

The DWT of these signals show that the gong and the dropped metallic object 
sounds have similar characteristics in the transform domain. This is due to their 
exponential type behavior in the time domain. However, the latch sound has a very 
particular transformation and the transition instances appear clearly throughout the scales. 
A preliminary analysis of the plots shows that the sounds produced by the dropped 
metallic object and the gong can be classified in the same group of transients, while the 
sound produced by the latch belongs to another group. It is possible to observe these 
features directly from the time domain plots, maybe due to a high SNR. Hence, a wavelet 
transform is not necessary. We have seen during the simulation phase that if the SNR is 
low, the transition instances become harder to identify from the time domain plots Using a 
DWT is a good technique to determine these transitions. 

Obtaining more metallic sound transients and transforming them, using the DWT 
for the purpose of their classification into groups, is recommended as a future study. 

2. Radar Pulses 

Four types of radar pulses were obtained from two types of radar at the Naval 
Postgraduate School. The first radar recorded was an AN/APS 31-A operating at a center 
frequency of 9.379 GHz, having a nominal pulse duration of 0.5 microseconds. The other 
three types of radar pulses were generated by an AN/SPS 64 (V) radar operating at a 
center frequency of 9.4 GHz and nominal pulse duration of 1.0, 0.5 and 0.15 
microseconds. All signals were heterodyned to 500 MHz using a local oscillator and a 
mixer, then digitized at 2 GHz with 8 bits per sample. There are 64 pulses per record, 8 
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records for the first type of radar pulse, 15 records for the second, 8 records for the third 
and 8 records for the fourth. Figures 36.a through 36.d show the first pulse of the first 
record of each one of the radars. Figures 37 through 44 show the last four scales of the 
DWT of six pulses of the first record of each radar using the Daub 12 wavelet. 

The DWT was applied to the first six pulses of each record. Although the pulses 
inside one record were generated sequentially by the same radar, their DWT scales do not 
present a particular repeating feature. Therefore, at this point, we could not classify the 
transformed radar pulses. Statistically, the number of radar pulses that were transformed 
was not enough to make any assessments. Hence, this is an area where additional research 
is recommended. 

3. Turn-On / Turn-Off of Push-To-Talk Radios 

These transient signals were generated by the turn-on and turn-off of six different 
models of transmitters produced by the same maker. They were collected and recorded by 
the Naval Security Group Activity, Charleston, SC. Ten samples of each of six 
push-to-talk Motorola radios were recorded. The carrier frequency of all of them was set 
to 138.525 MHz. The signals were passed through a 1 MHz bandwidth filter and then 
digitized with a sampling frequency of 5 MHz and a center frequency of 1.075 MHz 
Figures 45.a through 45.f show one sample of the tum-on of each transmitter and Figures 
46.a through 46.f show the corresponding transmitter turn-off .The last four scales of the 
DWT of the tum-on transients, using a Daub 12 wavelet, are shown in Fgures 47 and 48 
and the turn-off transients in Figures 49 and 50. 

The localization of the exact moment of the tum-on or turn-off of the transmitter 
can be extract from the DWT of these transients. Although these signals are sinusoidal, 
and as such the DWT should not perform well on them, their behavior in the transform 
domain is different before and after the switching instant. This might be due to the noise 
present in the records. Denoising the signals before applying the DWT is a possibility and 
deserves further studies. The transition instances can also be observed from the time 
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domain plots for most of the signals. If we observe Figure 45.d, the switching instant is 
difficult to identify, although if we look at Figure 48 the transition instance is observable at 
levels ten, 11 and 12. A spectrogram of the record in question was performed. Figures 
51. a and 51 .b show the result of the spectrogram using a hamming window and without a 
window, respectively. The instant of the switching time can also be observed in both 
spectrograms. The results of wavelet transforming these signals are inconclusive, since the 
spectrogram produces a good localization of the transition instant. A feature extraction 
algorithm based on the energy of each scale for the purpose of classification other than 
detection was implemented. The results were not satisfactory, since the levels of energy 
for a given scale seemed to be similar for all records and all transmitters. 
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Figure 40. DWT of Pulses 4, 5 and 6 of Radar AN/SPS-64 (V), 
Pulse Duration = 1.0 microsecond 
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Figure 42. DWT of Pulses 4, 5 and 6 of Radar AN/SPS-64 (V), 
Pulse Duration = 0.5 microseconds 


55 









































































































































































































































































































































































































































































3500 4000 


Figure 46. First Record of the Turn-Off of: 

(a) Transmitter 1 

(b) Transmitter 2 

(c) Transmitter 3 
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IV. CONCLUSIONS 


A. RESULTS 

The identification of transient signals using the Discrete Wavelet Transform 
(DWT) provided mixed results. When the signals of interest have sharp variations in the 
time domain and the Signal-to-Noise-Ratio is positive, it is possible to identify these 
behavior changes in the scales. The wavelet used as a basis function determines the overall 
performance of the algorithm. The use of the ad hoc wavelet as a basis function improved 
the identification of the test signals in low SNR, possibly because it matched the signal 
characteristics better and enhanced them relative to the noise. 

The DWT of the sounds produced by metallic objects indicated that it might be 
possible to use the wavelet transform for classification purposes. Also, the latch sound 
analysis showed that the exact instant of the transient occurrences can be established. 

The question of what should be the scale of interest was brought up during the 
simulation and the test phase. The higher scales have more energy and better time 
resolution and should; therefore, have more importance than the lower scales. 

A transient detection algorithm for the transmitters turn-on / turn-off based on the 
energy of the scales was implemented, but the results were inconclusive. The transformed 
signals have the same total energy at any given scale when compared between signals. This 
was perhaps due to the transmitters being manufactured by the same maker. 

The DWT is simple to implement and is very fast. Its speed depends on the number 
of scaling function coefficients that is being used. 

B. RECOMMENDATIONS FOR FUTURE STUDIES 

The implementation of an automatic transient detection algorithm based on the 
scales of the signal is the goal. Further studies, focused on finding an optimum wavelet 
function, are necessary. The use of a larger number of actual transient signals with known 
characteristics is recommended to make the conclusions statistically more reliable An 
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additional untapped source of information is the phase of the DWT, which should be 
investigated if the signals of interest are complex. 
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APPENDIX A. DEVELOPMENT OF NECESSARY CONDITIONS 


Necessary Condition 1: 

Integrating both sides of Equation 2.27 


J q>(t)dt = J 2 ho(k)(p(2t-k)dt 


= 2 h Q (k)\ cp(2 1 - k)dt. 


(A.l) 


Let z = 2t - k, so dz = dt/2, substituting into Equation A. 1 


J c p(t)dt =^S h 0 (k)j (p(z)dz , 


(A.2) 



(A3) 


The left hand side of Equation A.3 is equal to 1 and so 


^Z*o(*)=l. 

+ k 


or 


2 ho(k) = 2. (A.4) 

k 
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Necessary Condition 2: 

Substituting Equation 2.27 in Equation 2.39 we obtain 


<cp(0cp(r -k)) = J 2 ho(n)(p(2t - h 0 (m)(p(2(t -k)- m)dt 

n m 


= 2 2 /jo(h)/JoOOJ Cp(2r- «)(p(2r -2k- m)dt = d(k) . (A.5) 

n m 


The integral of Equation 2.44 is non-zero only if 

2t-n = 2t-2k-m => n - 2k + m. 

That leads to 

Jcp(2r-«)(p(2/ -2k-m)dt = -m-2k) . (A. 6) 

Substituting in equation A.5 we obtain 

2 2 Z»o(«)^o( w )i8(/i -m-2k) = 6(k) . (A.7) 

The Kronecker delta on the left hand side exists only if m-n-2k and Equation A.7 is 
now 

^I,ho(n)ho(n-2k) = m. 

or 

'Zh 0 (n)h 0 (n-2k) = 2h(k). (A.8) 

n 


70 




Necessary Condition 3: 
Reproducing Equation A. 8 


X ho(n)ho(n - 2k) = 2b(k ). 

n 

Let k = 0, and substitute in Equation A.8 

£ ha(n)h<i(n) = 28(0) = 2. (A.9) 

n 

Necessary Condition 4: 

Substituting Equation 2.47 and summing over £ we get 
X 2 ho(n)ho(n - 2k) = 2. (A. 10) 

k n 

Splitting Equation A. 10 in even and odd terms 

X X h 0 (2m)h 0 (2k + 2w) + X h 0 (2k + 2m + \ )h 0 (2m + 1) =2, (All) 

k L m 

rearranging terms 

X X^o(2A:+2/w) /to(2/w) + X X^o(2/r + 2/n+ 1) h 0 (2m+ 1) = 2 . (A 12) 

n? L ^ J #n L * 

Substituting Equations 2.41 and 2.42 into Equation A. 12 

SiX^o(2m) + S 2 X^o(2m+1) = 2. (A 13) 

m m 

Applying the same equations again in Equation A. 13 
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s\+sl = 2. 


(A. 14) 


Solving for S,and S 2 in Equations A. 10 and A. 13 we obtain 


Si=2>o(2«)=1. 

n 

Si = Yi hoQn + 1 ) = 1 . 

n 


(A. 15) 


(A. 16) 
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DWT of OOK, SNR = 3 dB 
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Figure B40. DWT of QPSK, SNR = 3 dB 
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APPENDIX C. MATLAB CODES 


mapdn.m 

Purpose: 

Computes a two-dimensional array A which defines the mean square map for a 
dimensional function/. 

Synopsis: 

A = mapdn(f,N,type). 

Description: 

Uses a = wavedn(f,N,type) to compute the wavelet transform of/ which is a sequence of 
length 2n; then reorders and arranges the square of elements ( l:2n ) to form the 
two-dimensional array A, which has n+1 rows and 2(n-l) columns; the volume under this 
surface is proportional to the mean-square value of f, using four mother wavelets: Haar, 
Daubechies, Spline and the ad hoc. 

Algorithm: 

The elements of the wavelet transform a are recorded according to the strategy described 
in Ref. 12, Chapter 17, so that when the matrix A is formed each element (squared) is 
located at the center of the wavelet it represents. 

Limitations: 

Restricted to even values of N in the range 2 to 20 for Daubechies’ wavelet functions, 
being N = 2 the call for the Haar wavelet function and N = 4 to the ad hoc wavelet 
function. No limitation for the number of spline function coefficients. 
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Modifications: 

Addition of two other types of wavelet functions to be used: spline and ad hoc. 
Plotting subroutine appended to the end of the program. 

function A = mapdn(f,N,type) 

% 

% Copyright (c) D. E. Newland 1992 
% All rights reserved 
% 

% Modified by Jorge Pitta, by permission from the author. 

% 

M = length(f); 
n = round(log(M)/log(2)) 
a = wavedn(f,N,type); 
b(l) = a(l);b(2) = a(2); 
for j = l:n-l 
for k= 1:2^ 
index = 2 / j+k+N/2-l; 

while index > 2 A (j+l),index = index-2^j;end 
b(index) = a(2^j+k); 
end 
end 
a-b; 

forj = 1 : 2 * 01 - 1 ) 

A(I j) -a(l); 
end 

forj = 2:n+l 
for k = 1:2 A (j-2) 
for m = l^n-j+l) 
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A(j, (k-1) *2 A (n-j+1 )+m) = a(2 A G-2)+k); 
end 
end 
end 

A = A.*conj(A); 

%% Printing Subroutine 

subplot((n+1 )/2,2,1 ),plot(abs(f),'r l ),title('signal-D AUB-COEFFS'), 
ylabel('f(x)') ,xlabel('x'); grid 
hold on 

for h = 2:l:n+l 

subplot((n+l)/2,2,h),plot(A(h,:),'r'),grid,ylabel(['level', num2str(h-2)]) 
end 

hold off 
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wavedn 

Purpose: 

Wavelet Transform. 

Synopsis: 

a = wavedn(f,N,type). 

Description: 

Returns an array whose elements are the wavelet transform of the sequence of elements in 
f. The analysing wavelet has N even coefficients. 

Algorithm: 

It is the Mallat’s Pyramid Algorithm adapted to a circular form. 

Limitations: 

It is necessary for the sequence f to have 2n elements, where n is an integer number. Its 
transform a also has 2n elements. The wavelet’s coefficients are imported from the file 
dcoeffs.m ( for the Haar and Daubechies wavelets), dcoeffsl.m ( for the Spline wavelet ) 
and dcoeffs2.m (for the ad hoc wavelet). 

If the sequence/is complex, then the transform a is also complex. 

Modifications: 

Addition of a third local variable called type that determines the type of wavelet function. 
Addition of three IF statements to use with the local variable type. 
function a = wavedn(f,N,type) 

% 

% Copyright (c) D. E. Newland 1992 
% All rights reserved 


124 



^ Modified by Jorge Pitta, by permission from the author. 
% 

M = length(f); 
n = round(log(M)/log(2)); 
if type = = 1 ;c = dcoeffs(N);end 
if type = = 2;c = dcoeffs 1 (N);end 
if type = = 3;c = dcoeffs2(N);end 
clr = fliplr(c); 

for j = 1:2:N-1 , clr(j) = -clr(j); end 
a = f; 

for k = n:-l:l 
m = 2 A (k-l); 
x =[0];y = [0]; 
for i = 1 :m 
forj = 1:N 
kO) = 2*i-2+j ; 

while kQ) > 2*m, k(j) = k(j)-2*m ;end 
end 

z = a(k); 

[mr,nc] = size(z); 
if nc > 1 , z = z’; end 
x(i) = c*z; 
y(i) = clr*z; 
end 

x = x/2 ; y = y/2 ; 
a(l:m) = x; 
a(m+l:2*m) = y; 
end 
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dcoeffs.m 


Purpose: 

Generates the Haar wavelet coefficients for N = 2 and the Nth order Daubechies’ wavelet 
coefficients up to N = 20. 

Synopsis: 
c = dcoeffs(N). 

Description: 

It is called by a = waden(f,N,type) for providing the N Daubechies wavelet coefficients. 

Limitations: 

Restricted to even values of N in the range 2 to 20. 

% 

% Copyright (c) D. E. Newland 1992 
% 

% All rights reserved 
% 

function c = dcoeffs(N) 

% 

ifN = = 2 
C' [11]; 
end 


ifN = = 4 

c = [(l+sqrt(3))/4 (3+sqrt(3))/4 (3-sqrt(3))/4 (l-sqrt(3))/4]; 
end 
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if N = = 6 

q = sqrt(10);s = sqrt(5+2*q); 

c = [(l+q+s)/16 (5+q+3*s)/16 (5-q+s)/8 (5-q-s)/8 (5+q-3*s)/16 (l+q-s)/16J; 
end 


if N = = 8 

c = [.3258030428051 ,1.010945715092 .892200138246,-.039575026236, - 
.264507167369,.043616300475,.046503601071,-0.14986989330]; 

end 


ifN = = 10 

c = [.226418982583,.853943542705,1.024326944260,. 195766961347,-.342656715382, 
-.045601131884,. 109702658642,-.008826800109,-.017791870102, .004717427938]; 
end 


ifN = =12 

c = [.157742432003, .699503814075,1.062263759882, .445831322930, 
-.319986598891,-. 183518064060,. 137888092974,.038923209708, 
-.044663748331 ,.000783251152,.006756062363,-.001523533805]; 
end 


ifN = = 14 

c = [.110099430746,-560791283626,1.031148491636,-664372482211, 
-.203513822463,-.316835011281,.! 00846465010,. 114003445160, 
-.053782452590;-.023439941565,.017749792379,.000607514996, 
-.002547904718,.000500226853]; 
end 


127 




ifN = =16 

c = [.076955622108,.442467247152,.955486150427,.827816532422, 
-.022385735333,-.401658632782,-000668194093,.182076356847, 
-.024563901046,-.062350206651,-019772159296,.012368844819,- 
.006887719256,-.000554004548,.000955229711 ,-.000166137261]; 
end 


ifN = 18 

c = [.053850349589,.344834303815,-855349064359,.929545714366, 

. 188369549506,-.414751761802,-. 136953549025,.210068342279, 
.043452675461 ,-.095647264120,.000354892813,.031624165853, 
-.006679620227,-.006054960574,.002612967280,.000325814672, 
-.0003 56329759,-000055645514]; 
end 


ifN = = 20 

c=[.037717157593,-266122182794,.745575071487,.973628110734,39763774177;. 
-.353336201794,.277109878720,. 180127448534,. 131602987102,. 100966571196, 
.041659248088,-046969814097,.005100436968,.015179002335,.001973325365, 
.002817686590,-.000969947840,-.000164709006,.000132354366,-.000018758416]; 
end 
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dcoeffsl.m 


Purpose: 

Generates spline wavelet coefficients for any integer N > 1. 

Synopsis: 
c = dcoeffsl(N). 

Description: 

It is called by a = waden(f,N,type) for providing the N spline wavelet coefficients. 

Limitations: 

None. 

% 

% Created by Jorge Pitta, 1994 

function c = dcoeffsl(N) 

% 

c = triang(N)'; 
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dcoeffs2.m 


Purpose: 

Generates the ad hoc wavelet coefficients for any integer N = 4. 

Synopsis: 
c = dcoeffs2(N). 

Description: 

It is called by a = waden(f,N,type) for providing the fourth order ad hoc wavelet 
coefficients. 

Limitations: 

Limited to only a fourth order ad hoc wavelet function. 

% 

% Created by Jorge Pitta, 1994 

function c = dcoeffs2(N) 

% 

c = [.5 1 .5 1]; 
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filtspec.m 


Purpose: 

Generates the spectrum response plot of the Haar, n* order Daubechies', spline and ad 
hoc wavelet function and their dilations. 

Synopsis: 

filtspec. 

Description: 

It calls subroutines tl.m ( for the Haar and Daubechies' wavelets ), t2.m ( for Spline 
wavelets), t3.m ( for the Ad Hoc wavelet). 

Created by Jorge Pitta, 1994. 

F = menu('WaveletVDaubechies','Spline','Ad Hoc','Own'); 
ifF = = l;tl;end 
if F = = 2;t2;end 
ifF = = 3;t3;end 


x=[l zeros(l,1023)];n=0:511; 

N=length(G); 
for nn=l:N 

H(nn)=((-1) A nn)* G(N+1 -nn); 
end 

w=filter(H, 1 ,x);e=filter(G, 1 ,x); 

W=abs(fftshift(fft(w, 1024)));E=abs(fftshift(fIt(e, 1024))); 
W1 =W(513:1024) ;E 1 =E(513:1024); 

H2=upsam(H,2); 
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Al=fft(H2,1024); 

Bl=fft(G,1024); 

C1=A1.*B1; 

H2G=fftshift(Cl); 

H2Gl=abs(H2G(513:1024)); 

H4=upsam(H,3); 

G2=upsam(G,2); 

A2=fft(H4,1024); 

B2=fft(G2,1024); 

C2=A2.*B2.*B1; 

H4G2G=fftshift(C2); 

H4G2G1 =abs(H4G2G(513:1024)); 

H8=upsam(H,5); 

G4=upsam(G,3); 

A3=ffi(H8,1024); 

B3=fft(G4,1024); 

C3=A3.*B3.*B2.*B1; 

H8G4G2G=fflshift(C3); 

H 8 G4G2 G1 =abs(H 8G4G2 G(513:1024)); 

H16=upsam(H,7); 

G8=upsam(G,5); 

A4=ffi(H16,1024); 

B4=fft(G8,1024); 

C4=A4.*B4.*B3.*B2.*B1; 

H16G8G4G2G=fftshift(C4); 

H16G8G4G2G1 =abs(H 16G8G4G2G(513:1024)); 
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H32=upsam(H,9); 

G16=upsam(G,7); 

A5=ffi(H32,1024); 

B5=fft(G16,1024); 

C5=A5.*B5.*B4.*B3.*B2.*B1; 

H32G16G8G4G2G=ffishift(C5); 

H32G16G8G4G2G1 =abs(H32G 16G8G4G2G(513:1024)); 

H64=upsam(H,l 1); 

G32=upsam(G,9); 

A6=fift(H64,1024); 

B6=ffi(G32,1024); 

C6=A6.*B6.*B5.*B4.*B3.*B2.*B1; 

H64G32G16G8G4G2G=ffishift(C6); 

H64G32G16G8G4G2Gl=abs(H64G32Gl 6G8G4G2G(513:1024)); 
m=max(H64G32G 16G8G4G2G1); 

y=linspace(0,pi,512);subplot(211) 

plot(y,H64G32Gl 6G8G4G2Gl/m,'k'),hold on 

plot(y,H32Gl 6G8G4G2Gl/m,'k'), 

plot(y,H 16G8G4G2G1 /m,'k'), 

plot(y,H8G4G2Gl/m,'k’), 

plot(y ,H4G2G 1 /m,'k'); 

plot(y ,H2G 1 /m,'k’); 

plot(y,Wl/m,'k') 

plot(y,El/m,’k’); 

hold off 
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tl 


Purpose: 

Calls the m.file dcoeffs.m for supplying filtspec.m with the Haar and Daubechies wavelet 
coefficients. 

Synopsis: 

tl. 

Description: 

Produces a pop-up menu on screen with the order of the Daubechies' wavelet to be used, 
if the order equals two, it calls Haar wavelet coefficients. 

Created by Jorge Pitta, 1994. 


O = menu('Order7274767871 0 ’,' 12 ’,' 1 4',’ 1 6',' 1 8',' 20 '); 

ZZ = 2*0; 

G = dcoeffs(ZZ)/sqrt(2); 
for zz = 1. 7.7 

H(zz) = ((-1 ) A zz)* G(ZZ+1 -zz); 
end 

KK = 'Daub'; 


t2 


Purpose: 

Calls the m.file dcoeffsl.m for supplying filtspec.m with the Spline wavelet coefficients. 
Synopsis: 

a. 


Description: 

Produces a pop-up menu on screen with the order of the Spline wavelet to be used, if the 
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Created by Jorge Pitta, 1994.. 


O = menu( , OrderV2V4V6V8','10 , , , 12V14 , , , 16V18 , , , 20 , ); 
ZZ = 2*0; 

G = dcoeffs 1 (ZZ)'/sqrt(2); 
for zz = 1 :ZZ 

H(zz) = ((-1 ) A zz)* G(ZZ+1 -zz); 
end 

KK = ’Spline'; 

t3 


Purpose: 

Calls the m.file dcoefFs2.m for supplying filtspec.m with ad hoc wavelet coefficients. 


Synopsis: 

t3. 


Description: 

Produces a pop-up menu on screen with the order of the ad hoc wavelet to be used. 
However, it is limited to an order 4. 


Created by Jorge Pitta, 1994. 
ZZ = 4; 

G = dcoeffs2(ZZ)/sqrt(2); 
for zz = 1 :ZZ 

H(zz) = ((-1 ) A zz)*G(ZZ+l -zz); 
end 

KK = ’Weird'; 
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